function [y1 y2] = funcKxS(nx, nz, zeta, CONSTS)

    k0 = CONSTS.k0;
    c = CONSTS.c;
    Z0 = 4*pi/c;
    eps = CONSTS.eps;
    eta = CONSTS.eta;
    
    y1 = (1i*Z0*k0/(2*pi))*(1+(2*nx.^2)./(1+abs(eps*eta))).*...
         ((sqrt(nz.^2+nx.^2).*cos(k0*nz.*abs(zeta)))./(nz.^2));
     
    y2 = (1i*Z0*k0/(2*pi))*(-2i*nx.^2./(sqrt(abs(eps*eta)))).*...
         ((cos(k0*nz.*abs(zeta)))./(sqrt(nz.^2-abs(eps/eta)*nx.^2)));
     
    if(false)
        figure; plot(nz, real(y1), 'b-', nz, real(y2), 'r-');
    end
    
end